Submitted to ApJ 



Resonant inclination excitation of migrating giant planets 

Edward W. Thommes 
Astronomy Department, University of California, Berkeley, CA 94720 

and 

Jack J. Lissauer 

Space Science Division, MS 245-3, NASA-Ames Research Center, Moffett Field, CA 94035 

ABSTRACT 

The observed orbits of extrasolar planets suggest that many giant planets mi- 
grate a considerable distance towards their parent star as a result of interactions 
with the protoplanetary disk, and that some of these planets become trapped in 
eccentricity-exciting mean motion resonances with one another during this mi- 
gration. Using three-dimensional numerical simulations, we find that as long as 
the timescale for damping of the planets' eccentricities by the disk is close to or 
longer than the disk-induced migration timescale, and the outer planet is more 
than half the mass of the inner, resonant inclination excitation will also occur. 
Neither the addition of a (simple, fixed) disk potential, nor the introduction of 
a massive inner planetary system, inhibit entry into the inclination resonance. 
Therefore, such a mechanism may not be uncommon in the early evolution of a 
planetary system, and a significant fraction of exoplanetary systems may turn 
out to be non-coplanar. 



1. Introduction 

The orbital migration of bodies locked in mean-motion resonances is a major recurring 
theme in the early history of the Solar System. A review of resonance dynamics in the 
Solar System is given by Peale (1976). Goldreich (1965) showed that tidally-induced orbital 
migration is likely responsible for most of the mean-motion commensur abilities among the 
satellites of Jupiter and Saturn. In the Uranian satellite system, tidally-induced convergent 
migration of the moons Miranda and Umbriel could have led to capture into, and subsequent 
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escape from, inclination-type resonances, thus producing Miranda's anomalously high orbital 
inclination of 4.3° (Tittemore & Wisdom 1989; Malhotra & Dermott 1990). 

In the case of the outer Solar System, Malhotra (1998) and Gomes (2000) showed that 
bodies entrained in the 2:3 exterior mean-motion resonance of an outward- migrating Neptune 
can also encounter the Kozai resonance, as well as secular and mean-motion inclination 
resonances. These resonances increase inclination, and this may be the mechanism by which 
Phito and the Kuipcr belt objects presently sharing the 2:3 resonance (plutinos) attained 
their high inclinations. 

Snellgrove et al. (2001) and Lee & Peale (2002) performed in-depth investigations of 
resonant capture and migration in the context of the GJ 876 planetary system. They found 
that convergent orbital migration of planets typically leads to their capture into the 2:1 mean- 
motion resonance. The resonant lock thus produced is characterized by the libration of the 
first-order resonance angles (cosine arguments in the expansion of the disturbing function; 
see for example Murray & Dermott (1999)), 



about zero degrees, where Aj is the ith planet's longitude, uji is its longitude of pericenter, 
and i = 1,2 refer to the inner and outer planet, respectively. Since Oe^ — Oe^ = ui — U2, the 
zeroth-order secular resonance angle 



also hbrates about zero; in other words, the orbits of the two planets librate about apsidal 
alignment. 

The analysis of Lee & Peale (2002) is restricted to the planar case, thus the effects 
of migration in resonance on orbital inclinations are not considered. Starting with slightly 
non-coplanar orbits, we initially see the same behaviour as reported by Lee & Peale (2002): 
In the absence of any damping mechanism, the eccentricities of both bodies grow. However, 
we find that once the inner planet's eccentricity is high enough, the system will also enter an 
inclination-type resonance, which induces rapid growth in the inclinations of both planets. 

In this article, we report the results of three-dimensional simulations of migrating planets 
in resonance. We give a summary of the underlying physical picture in Section 2. The 
methods that we use, as well as initial conditions, are presented in Section 3. Baseline 
cases with two equal mass planets are discussed in Section 4. We extend our analysis to 
planets of unequal masses in Section 6, impose eccentricity damping in Section 7, include 



= Al — 2A2 + Oi and 
Oe^ = Al - 2A2 + 0)2 



(1) 



6cj = 0Ji- U2 



(2) 
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the gravitational potential of the disk in Section 8, and add four Uranus/Neptune mass 
inner planets in Section 9. We conclude in Section 10 with a summary of our results and a 
discussion of their implications. 

2. The physical picture 

Gas giant planets must form in a protoplanetary disk which has not yet lost its gas; this 
constrains the time of formation to within, approximately, the first 10^ years of the disk's 
lifetime (Strom, Edwards, & Skrutskie 1990; Haisch, Lada, & Lada 2001). If gas giants form 
by the accretion of a massive atmosphere onto a ~ 10 solid core (Mizuno, Nakazawa, 
& Hayashi 1978; Pollack ct al. 1996; Bodcnhcimer et al. 2000), then the most favourable 
location for such large cores to form by planctcsimal accretion is in an annulus several AU 
wide beyond the ice-condensation line, roughly the location of Jupiter and Saturn in our own 
Solar System (Thommes, Duncan, & Levison 2003). A gas giant-sized body (mass of order 
10^ Mq) will open a gap in the gas disk. Its orbital evolution will subsequently be locked 
in to the viscous evolution of the gas disk, as long as the disk mass is large compared to 
that of the planet. This is commonly referred to as Type II migration (Ward 1997); see Lin 
et al. (2000) or Thommes & Lissauer (2003) for a review. In this picture, significant Type 
II migration probably did not occur in our Solar System; however, this mechanism is likely 
to have played a significant role for those extrasolar giant planets with orbital semimajor 
axes of less than a few AU. For a review of discovered extrasolar planets and the statistics 
of their orbits, see Gumming et al. (2003). 

If a disk forms two gas giants, the annulus of gas between their individual gaps may be 
relatively rapidly depleted, so that both planets end up sharing a common gap (Kley 2000; 
Bryden et al. 2000) Also, if the viscosity is higher and thus the accretion of the disk onto 
the star faster at smaller radii, as in the model of Gammie (1996), the inner disk may be 
depleted on a timescale short compared to the migration timescale of the planets. 

3. Numerical methods 

We use a variant of the SyMBA integrator (Duncan, Levison, & Lcc 1998), which has the 
desirable properties of the mixed- variable symplectic (MVS) mapping of Wisdom & Holman 
(1991), and in addition handles close encounters among objects by employing a variable 
timestep technique. We impose an inward migration on the outer planet, in a way which 
simulates interaction with a viscously evolving outer gas disk: We define an inner radius 
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'"disk to the gas disk; this radius is made to change at a rate rdisk- The semimajor axis of the 
outer planet, 02, is then made to change based on its location relative to the gap edge: 

. ^ r --S'(a2 - rdisk), 02 > rdisk 
^ I 0, 02 < rdisk 

where S is large and positive, so that the magnitude of 0,2 increases steeply outside of rdisk- 
In this way, a planet with 02 > rdisk is always constrained to move inward at a mean rate 
rdisk; if the torque required to maintain this migration rate changes (as happens when the 
planet resonantly captures an interior body which it then has to push along), the planet's 
position relative to the gap edge simply re-adjusts until its orbital decay rate again matches 
fdisk- It should be noted that, strictly speaking, rdisk corresponds not to the edge of the 
disk itself but rather to the location outwards of which a planet is subject to a nonzero disk 
torque. 

The speed with which the disk edge moves inward corresponds to the radial accretion 
speed of the disk. This is given by 

racer ~ ClQ;(rQ)(c/rQ)^ (4) 

where Q, is the Keplerian angular frequency, c is the gas sound speed, Ci is a constant of 
order unity, and a is the Shakura-Sunyaev viscosity parameter: 1/ ~ ac^/fl (e.g.. Ward 1997). 
The factor (c/rQ) is also the aspect ratio h of the disk, where h = H/r; H is the disk scale 
height at radius r. For simplicity we adopt the (flared) aspect ratio of the Hayashi (1981) 
nebula model, which has h oc r^/^-^ Chiang & Goldreich (1997) obtain a similar r-dependence, 
h oc r^/'^, for r ;^ 80 AU. We thus have a disk accretion speed racer which does not vary with 
r, assuming a is constant. Using h = 0.05(r/AU)^/'^, the accretion speed is racer ~ lO^^a 
AU/year. The a viscosity of a protoplanetary disk is estimated to be 10^"^ to 10~^ (^-g-, 
Cabot et al. 1987; Dubrulle 1993). Adopting a = 10^^, we have racer ~ 10^^ AU/year. 

It should be emphasized that migration is imposed on the outer planet as per Eq. 3 
regardless of the eccentricity of the planet, or its inclination relative to the plane of the disk 
(which is the planet's original orbital plane). Hydrodynamic simulations of the interaction 
of a signiflcantly eccentric and/or inclined gap-opening planet with a disk have, to our 
knowledge, not yet been performed, so we operate under the simple assumption that, no 
matter how the planet's orbit evolves, its semimajor axis remains locked relative to the inner 
edge of the disk. However, at sufficiently high eccentricity and/or inclination, a planet will 
very likely no longer be an effective barrier to the disk, so that this assumption would no 
longer be valid. We will discuss implications in Section 10. 

The shortest base timestep used in the simulations is 10~^ years, and integration is 
stopped when a body comes within 0.1 AU of the star. Beyond this point, the timestep is 
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deemed too long, and thus the azimuthal motion of the inner planet during one timestep too 
large, to accurately resolve the interactions between the two planets. 

In all of the simulations below, we start with a system of two planets on nearly circular 
orbits about a one solar mass (M©) star, the inner at 5 AU, and the outer at 8.5 AU. Initial 
eccentricities and inchnations are of order 10~^ and (10~^)°, respectively. 

4. The baseline case 

We begin with a system consisting of two Jupiter-mass planets. Several runs are per- 
formed, with differing initial orbital phases, but the outcomes are very similar. Fig. 1 shows 
one example. The outer planet migrates inward until it and the inner planet are mutu- 
ally captured into the 2:1 mean-motion resonance, with both first-order eccentricity- type 
resonance angles, 0^^ and (Eq- l); librating about zero, which means that conjunc- 
tions are occurring when both planets are at pericenter, and thus that the orbits are also 
librating about apsidal alignment. This behaviour is to be expected even for planets that 
start with a larger separation in orbital periods; Lee & Peale (2002) find that as long as 
a'2 < 10^'^(a/AU)^^/^, capture into the 2:1 resonance occurs with certainty if the eccentrici- 
ties of both planets are low (< 10~^), whereas larger eccentricities are required for a nonzero 
capture probability into more distant, higher-order resonances which would be encountered 
first (such as 3:1). 

The eccentricities of both planets increase as they migrate in resonance, which can be 
understood in terms of energy {E) and angular momentum (L) exchange between the disk 
and the system of resonantly-locked planets (Lissauer, Peale, & Cuzzi 1984): In order to 
change the semimajor axis of a circular orbit while keeping it circular, E and L have to be 
changed in the right ratio: E/L = dE/dL = a/ GM^/a?, which is just the orbital angular 
velocity, n{a). But n{a) is a decreasing function of semimajor axis. Thus, even though we 
assume a disk torque that changes the orbit of the initially circular outer planet in such a 
way that E/ L = ^(02), this ratio is too low for the inner planet. Since E and L are negative, 
the two-planet system as a whole acquires an angular momentum deficit and must respond 
by increasing the eccentricity of one or both of its components. The rate of eccentricity 
growth can be reduced by making the applied E/L higher than what would be required to 
keep the outer planet by itself circular. This would amount to applying eccentricity damping 
to the outer planet. The evolution of the system with eccentricity damping by the disk is 
investigated in Section 7. 

In Fig. IE and F, we also plot the resonance angles for the 4:2 inchnation-type mean 
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motion resonance. This is a second-order resonance (which is the lowest possible order for 
an inclination resonance); the associated terms in the expansion of the disturbing function 
are of order P. These angles are defines as 

e.2 = 2Ai - 4A2 + 2^1 (5) 

and 

Oq = 2Ai - 4A2 + 2^2, (6) 

with Qi, ^2 being the longitudes of the ascending nodes of the inner and outer planet, 
measured relative to an arbitrary but fixed direction in the inertial frame, and in the plane 
of the (notional) protoplanetary disk, which is, within an inclination of order (10~^)°, also 
the initial plane of both planets' orbits. 

Both angles are initially circulating, but ^,^2/2 and decrease over time until, at 

about 4 X 10^ years, 6q j1 starts librating about about 90°, while simultaneously Qi^jl starts 
librating about 270°. This means that ftj2 and Q^i^ are both librating about 180°. Fig. 
ID shows that initially, the inclinations of both bodies remain at < 0.01°, oscillating but 
undergoing no net growth. However, as the system enters the inclination resonances, the 
inclinations of both bodies begin to grow rapidly to ~ 10° (the outer planet) and ~ 20° (the 
inner one) in just tens of thousands of years. Thereafter, the inclinations continue to grow 
more slowly. 

This simultaneous hbration of and Q^^ also imphes the libration of the "mixed" 
resonance angle, 

Q,,,, = 2Ai - 4A2 + ^1 + ^2 = Qiil2 + ^,2/2 (7) 
about 0°, as well as hbration of the zeroth-order secular resonance angle, 

^^, = ^1-^2 = ^^2/2-^,2/2 (8) 

about 180°. The latter means that the two orbits are librating about anti-alignment of their 
fines of nodes. Fig. IH shows that the hbration ampfitude is only a few degrees. 

Furthermore, simultaneous libration of 6'ei, 6*62, and also implies the libration of the 
resonance angles 

^e,iii2 = Ai - 2A2 + cDi ± ((^1 - ^2) = 0e^ ± On (9) 

and 

9e,i,i, = Xi- 2A2 + CD2 ± (fii - ^2) = 0e, ± On. (10) 

The near-perfect anti-alignment of the lines of nodes can be understood in terms of the 
conservation of that component of the system's angular momentum which lies in the original 
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Fig. 1. — Evolution of the two planets in the baseline run, in which an inward migration rate 
of d = 10~^ AU/year is imposed on the outer one. Both planets are one Jupiter mass (1 Mj); 
quantities related to the inner (Planet 1) are plotted in blue, and in red for the outer (Planet 
2). The panels show: A. Semimajor axis, periastron distance and apastron distance; B. The 
lowest-order eccentricity-type resonance angles (blue) and (red); C. Eccentricities of 
both planets; D. Inclinations of both planets, as well as (black) the mutual inclination; E., F. 
6*^2/2 and 0^2/2 where 9^2 and are the lowest-order inclination-type resonance angles; G. 
The angle between the planets' longitudes of periastron; H. The angle between the planets' 
lines of nodes. 
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Fig. 2. — Evolution of the planets when the migration rate imposed on the outer planet is 
10~^ AU/year, one tenth the rate in the baseline run (Fig. 1). All panels are as in Fig. 1. 
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orbital plane. Since this is initially very near to zero, the in-plane components of the two 
orbits' angular momentum vectors have to be pointing in very nearly opposite directions 
once the inclinations get large, otherwise they cannot sum to (almost) zero. This also means 
that the relative inclination of the two planets' orbital planes (see ID) is the sum of their 
individual inclinations. Thus the relative inclination of the two planets very quickly goes up 
to 30°. 

At 5.2 X 10^ years, the system evolves out of both inclination-type resonances. There 
is an abrupt increase in inclinations at this time, but thereafter, as is to be expected, there 
is no further net inclination increase. At the same time, the inner planet's eccentricity-type 
resonance angle, Oe^, switches from libration about 0° to "near-libration" about 180°: The 
angle librates about 180° for a few cycles, then makes a complete revolution, librates for 
another few cycles, and so on. Prom a perturbative analysis point of view, this means that 
the contribution from the corresponding term in the expansion of the disturbing function 
continues to be important; its cosine argument spends most of its time around 180° so the 
term does not average to zero over time. An interesting point is that in a completely planar 
system, this transition in libration center is inevitably cut short by the system becoming 
unstable. This happens because the system is trying to evolve from conjunctions at both 
planets' pericenters, to conjunctions at the outer planet's pericenter and the inner planet's 
apocenter. Somewhere in between, then, the two bodies will reach a conjunction at the point 
where their orbits cross, bringing about a very close encounter or even a physical collision 
(Lee & Peale 2003, and Lee & Peale, private communication). However, if a significant 
inclination between the two orbits has been attained prior to this transition, as is clearly the 
case in Fig. 1, this will serve as a protection mechanism against close encounters, allowing 
the system to stably evolve through the transition in ^ei's libration center. 

The simulation is stopped at a httle after 7x10^ years, when the inner body comes within 
0.1 AU of the star. Other runs with randomly varied initial orbital phases, as well as different 
intial positions of the planets (though still starting outside the 2:1 mean-motion resonance) 
produce a similar evolution; in each case the the inclination resonance is encountered after 
capture into the 2:1 eccentricity resonance, and inclinations grow rapidly thereafter. 

5. Dependence on migration rate and initial inclination 

Another example is shown in Pig. 2. In this case, we reduce the speed at which the 
disk edge moves inward, racci , by an order of magnitude, which amounts to letting a = 10~^. 
The evolution continues to be very similar to that shown in Pig. 1, scaled up by a factor of 
ten in time. However, in slowing the migration rate we have also made the migration more 
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adiabatic; that is, we have reduced the amount of migration during one hbration period 
of the resonance(s). Thus, it is to be expected that the dynamics change by more than a 
simple rescahng in time. One difference apparent from Fig. 2 is that hbration amphtudes 
are reduced; the average values of e and i, however, are very close to those in the faster 
migration. Another difference is that the times at which the inclination resonances are 
entered and exited do not scale exactly with the slower migration speed. 

We also probe the dependence of the system's evolution on the initial mutual inclination 
of the two orbits. Up to now, we have made this ~ 0.01°. We perform a run where the mutual 
inclination is decreased to ~ 0.0005°, and another where it is increased to 1°. The resulting 
inclination evolution of the two planets is plotted in Fig. 3, together with the baseline case 
of the previous section. The smaller initial inclination causes very little change in the time of 
onset of the inclination excitation, the timescale over which the inclination is subsequently 
excited (~ 2 x 10^ years), or the value to which the mutual inclination is subsequently excited 
(~ 30°). The larger initial inclination causes a slightly earlier onset of inclination growth, 
by an amount At/t ~ 0.05. Evidently, therefore, entry into the inclination resonance is not 
strongly influenced by the actual value of the initial mutual inclination, as long as it is small. 



6. Varying the mciss ratio 

We next investigate the effect of making the mass ratio of the two planets, M1/M2, differ 
from unity. We find that when M1/M2 ^ 2, the inclination resonance is no longer reached. 
This is because onset of the resonance requires the inner body to have a large eccentricity; 
however, when there is a large disparity in the masses of the two planets, the eccentricity 
of the more massive one will not grow very much. Fig. 4 shows an example in which 
M1/M2 — 3. Although the outer body's eccentricity has reached 0.85 by the end of the run, 
that of the inner body does not get above 0.2, and the inclination resonance is avoided. We 
defer a detailed investigation of the conditions for onset of the inclination resonance to future 
work. For the time being, it is useful to recall that the simultaneous hbration of a body in 
the 2:1 eccentricity-type resonance and the 4:2 inclination-type resonance requires that its 
longitude of pericenter precession rate, (D, be on average equal to the precession rate of its 
line of nodes, Ct (as in a Kozai resonance, e.g., Murray & Dermott 1999). Both of these rates 
change as the eccentricity increases; for a body already in the eccentricity-type resonance, 
we expect that the inclination-type resonance sets in when the rates become equal. 

When the inner body is smaller, the inclination resonance is still entered for much more 
disparate masses. Fig. 5 shows the evolution when the outer planet is Jupiter-mass and the 




Fig. 3. — Evolution of the planets' inclinations for differing initial mutual inclinations of 
< O.OOr (top), O.Or (middle; the baseline case), and 1° (bottom) 
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inner is Earth-mass, making M1/M2 ~ 1/300. Though the behaviour is quahtatively similar 
to the equal-mass case (Fig. 1), there are some obvious differences: As is to be expected from 
conservation of momentum, the small planet suffers far more inclination excitation than the 
large one. Also, there is a delay of ~ 5 x 10^ years between the onset of libration in 6^2 and 
9^2 , and the inclination growth, which begins with the libration of , is not as fast as in the 
equal-mass case. 

We perform additional simulations in which M1/M2 ~ 10^^, making it effectively a test 
particle (with M2 still being Jupiter mass. Mi is now comparable to the largest asteroids). 
Yu & Tremaine (2001) find that a test particle in a 2:1 resonance with an inward-migrating 
massive body, initially inclined relative to the massive body by of order 1°, can have its 
inclination excited all the way up to 90° and then be released from the 2:1 resonance alto- 
gether. When we use a similar initial inclination, we get capture into inclination resonance 
and excitation of the small body's inchnation. We only follow the evolution of the system 
until the inner body's pcriccntcr drops below 0.1 AU; at this point the body is still trapped 
in both the eccentricity and inclination resonance, and its inclination, which is still growing, 
has reached ~ 70°. It is now only and Oq, not and Oq, which librate. With an initial 
mutual inclination of 0.01°, however, it is possible for the small body to quickly pass through 
the inclination resonance, "librating" for less than a full cycle before leaving it again, while 
having its inchnation excited to only ~ 1°. Though Yu & Tremaine (2001) do not check 
for libration in inclination- type resonances, the similarity in evolution of eccentricity and 
inclination between our simulations and their Fig. 5 leads us to believe that the behaviour 
they see is the M1/M2 — >■ hmit of the mechanism we explore here. In any case, we do not 
investigate the behaviour at such extreme mass ratios in more detail in the present work; if 
we are to consider pairs of giant planets, the possible masses range from 10^ to 10"^ M®, so 
the smallest possible ratio is of order M1/M2 ~ 10~^, as in the case pictured in Fig. 5. 



7. Damping of eccentricities 

Next, we consider the possibility that planet-disk interactions may damp planetary 
eccentricities. Following Lee & Peale (2002), we make the relative eccentricity damping rate 
proportional to the relative applied semimajor axis decay rate for the outer planet: 

^=5^. (11) 
62 0-2 

The applied decay rate, a^^^, is the rate at which the outer planet by itself would migrate 
under the action of the applied torque. However, because the outer planet has to push along 
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Fig. 4. — Evolution of the planets in the case where the inner is Jupiter-mass, and the outer 
is one third that, i.e., Saturn-mass (Section 6). The panels are as in Fig. 1. 
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Fig. 5. — Evolution of the planets in the case where the outer planet is Jupiter-mass, and 
the inner ~ 1/300 Jupiter mass, i.e.. Earth-mass (Section 6). The panels are as in Fig. 1. 
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the resonantly-locked inner planet, its actual migration rate, d2, is less than . Recall 
that torque is applied in such a way as to keep d2 (nearly) constant (Eq. 3). 

This eccentricity damping rate is applied to the outer planet only; the inner is, again, 
assumed to not be interacting with the disk directly. We perform simulations with S — 2 
and S — b. Results are shown in Fig. 6. In the S — 2 case, the inchnation resonance is 
entered later due to the slower eccentricity growth, and the maximum mutual inclination 
reached is less than 20°. When S = 5, the planets' eccentricities reach equilibrium values 
before the inclination resonance is attained, and the system thus remains coplanar. 

One effect we have not modeled is direct damping of the planets' inclinations by the 
disk. Lubow & Ogilvie (2001) find that such damping should occur when a > 10~^, and 
proceed on a timescale of order 10^ years. Thus, our baseline runs, with a = 10~^, are at 
the limit of where damping is expected to occur. Also, the inclination damping timescale 
is long compared to the timescale of inclination growth when the inclination resonances are 
first encountered, which is only a few x 10^ years. Thus, damping of inchnations by the gas 
disk likely has little effect on inclinations excited by resonances among planets. 



8. The disk potential 

Thus far we have only included the effect of the gas disk in the form of a forced migration 
experienced by the outer planet. However, the disk may inffuence the dynamics of the 
planets in other important ways. Specifically, the gravitational force from a disk will modify 
the planets' nodal and apsidal precession rates, Cl and a). This can change the locations 
and relative spacings of the mean-motion resonances, and so might change the nature of 
the system's resonant evolution. Therefore, it is important to incorporate the potential of 
the disk into our analysis. We assume that the "disk" is a thin annulus of surface density 
S(r) located in the z = plane, with inner radius ri and outer radius r2, where r is the 
stellocentric distance. The potential due to the disk, as felt at a radius R and height z, is 

*^^'") = "/ / /o2^ 2! 2 9 ^^' ^12) 
^0 v it"' H 2;^ — 2rR cos 

and the corresponding force is F^iR^z) = —d^/dR, Fz{R, z) = —d^/dz. To obtain the 
force due to an outer disk on a body, we can expand the numerator in powers oi R/r and 
z/r, then perform the integration. Assuming a surface density of the form 

E(r) = Eor-^/^ (13) 



i.e., having the density profile of the standard Hayashi (1981) nebula model, we obtain 
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Fig. 6. — Evolution of the two-planet system when the eccentricity of the outer is damped. 
The left set of panels show the case where the relative damping rate, 62/62, is twice the applied 
relative semimajor axis decay rate, a^"^^ /a2'-, the right set show the case where the relative 
damping rate is 5x a^^^/a2. The individual panels show: A. Semimajor axis, periastron 
and apastron distance; B. Eccentricities; C. Inclinations; D. The mixed inclination-type 
resonance angle Oi^i^ 
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and 



r=r2 

(14) 



r=ri 



(15) 



r=ri 



We performed a simulation with a disk surface density of 5 times E^j^, where 

= 1.7 X 103(r/l AU)-3/2 g/cm' (16) 

is the Hayashi minimum-mass surface density. The inner radius of the disk is made to be 20 
Hill radii of the outer planet beyond raisk at all times. The Hill (or Roche) radius of a planet 
of mass M orbiting at a distance r from a primary of mass M^, is th = (M/3M*)^/"^r. The 
actual width of a planet-induced gap is likely significantly smaller, between 3 and 4 Hill radii 
(e.g., Bryden, Rozyczka, Lin, & Bodenheimer 2000). However, since the migrating planets 
tend to acquire large eccentricities, we endeavour to keep the disk edge far enough away that 
it lies well beyond the planets' apocenters at all times, thus keeping the expansion of the 
potential valid. The response of a disk to an eccentric gap-opening planet is a complicated 
problem which has thus far not been investigated (e.g., Lin et al. 2000), so to avoid having to 
make arbitrary assumptions, we restrict ourselves to modeling the gravitational interaction 
of the planets with the more distant part of the disk. 

The evolution of the system is shown in Fig. 7. Comparing to Fig. 1, one can see that the 
evolution is quite similar to the diskless case. Despite the additional radially-varying apsidal 
precession induced by the disk potential, the two planets again become apsidally aligned 
upon encountering the 2:1 resonance. Also, despite the change in the vertical frequency, the 
inclination resonance sets in at nearly the same time as in Fig. 1. One noticeable difference is 
that the angle between the planets' fines of nodes, ^21 — ^2, has a larger amplitude of libration 
about 180°, and in fact circulates briefly several times. This is possible because, with the 
addition of the disk potential, the two-planet system no longer conserves the component of 
its angular momentum which lies in the disk plane. 



9. Additional planets 

Finally, we investigate the effect that additional smaller planets have on migration in 
resonance. We add four bodies, each 15 in mass (i.e., Uranian or "ice giant" size), 
interior to the inner of the two Jupiter-mass bodies. The bodies are separated by between 
10 and 15 mutual Hill radii, so as to form a relatively stable system. The innermost initially 
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Fig. 7. — Evolution of the two-planet system with the addition of the potential of a 5x 
minimum- mass outer disk. Panels are as in Fig. 1. 



-19- 



has a semimajor axis around 1 AU, while the outermost is at 3 AU, just short of the interior 
2:1 resonance of the Jovian body at 5 AU. We do not argue that it is hkely for a planetary 
system's inner region to contain so many large bodies at any time in its evolution. In-situ 
formation alone is, at any rate, unlikely to produce such massive bodies in this region (e.g., 
Lissauer 1987; Thommes, Duncan, & Lcvison 2003). Rather, our intent is to put into the 
simulation about as much mass in sub-gas giant bodies as will stably fit in the region. In 
this way, we obtain an idea of the upper limit (apart from stochastic close approaches) on 
the perturbations that the inner part of the planetary system could exert on the resonant 
migration of the two gas giants. 

Fig. 8 shows the resulting orbital evolution, including semimajor axes of the Uranian- 
mass bodies, in one of a set of eight simulations performed. The last of these extra bodies is 
eliminated at 3x10^ years. All are removed by crossing the inner simulation boundary, which 
in these runs is at 0.5 AU. Comparing to Fig. 1, there are some key differences between 
the orbital evolution of the two gas giants here and in the baseline simulation. Here, the 
onset of the inclination resonance, shown by the transition of 6q and 6j2 from circulation to 
libration, takes place about 2 x 10^ years later than in the baseline simulation. However, the 
gas giants already acquire significant inclinations before this happens. The growth of the 
inclinations begins at about the same time as 9^, the difference between the lines of nodes, 
begins to librate about 180°. This means that, although the two planets are not in the pure 
inclination-type resonances, they are (since and 6*62 are librating at the same time) in the 
mixed eccentricity-inclination resonances involving the libration of 9eini2 and 9e2hi2- 

The pure inclination-type resonances do not set in until after the last Uranian planet is 
eliminated, and final inclinations are comparable to those reached at the same time in the 
basehne run. The transition from hbration about 0° to libration about 180° of one of the 
eccentricity-type resonances seen in Fig. 1 does not occur here, simply because the inner 
boundary in this run is at 0.5 AU instead of 0.1 AU (to allow a longer timestep and thus 
speed the integration) , and therefore the inner Jovian planet is eliminated before its libration 
center can switch. 

With this many bodies, the system becomes strongly stochastic. Of the eight runs 
performed, six evolve in a manner qualitatively similar to the one shown in Fig. 8: The gas 
giants undergo resonant capture and evolve inward together, pushing the interior planets 
ahead of themselves. In all of these cases, significant resonant excitation of the giant planets' 
inclinations occurs, though in some cases a 3:1 rather than a 2:1 commensurability is involved. 
In the two remaining runs, the gas giants undergo a physical collision. This happens because 
close encounters with the smaller planets can exert a strong enough perturbation on the large 
planets to temporarily knock them out of resonance with each other. The large planets' orbits 
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Fig. 8. — An example of the evolution of the system with four 15 (Uranian mass) planets 
added interior to 3 AU. Panels are as in Fig. 1, and additionally the semimajor axes of the 
Uranian planets are shown (in black, green, cyan and magenta) in Panel A. The last of them 
is eliminated at about 3 x 10^ years. 
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Fig. 9. — The full set of eight runs performed with four inner Uranian planets, showing the 
evolution of the two Jupiter- mass planets' inclinations in each case. In runs 1 and 3, the two 
Jovian planets merge before the inclination resonances are encountered. 
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already cross by this time, so while their resonant lock is removed, they are free to undergo 
close encounters with each other. 

Fig. 9 shows the evolution of the inclinations in each of the eight runs. Runs 1 and 
3 are the ones where the gas giants collide, in both cases before significant inclinations are 
excited. In all of the other six runs, relative inclinations in excess of 30° arise. 

10. Discussion 

We have shown that two planets migrating in resonance can undergo significant exci- 
tation of their inclinations under a wide range of conditions. The excitation mechanism 
is libration in the lowest order (order i^) inclination-type mean-motion resonances. Even 
though convergent migration initially leads to capture into a first-order eccentricity-type 
resonance, at sufficiently large eccentricity, the two bodies can librate simultaneously in 
both types. Previous work has also suggested the possiblity of high inclinations arising dur- 
ing the formation of planetary systems (Levison, Lissauer, & Duncan 1998; Yu & Tremaine 
2001). In the former case, the inclination excitation mechanism considered was stochastic 
close encounters among planets. However, since these typically eventually lead to ejections 
and mergers, one would expect few long-lived non-coplanar systems of closely-spaced planets 
to be produced in this way. In the latter case, the authors showed that massless test particles 
resonantly trapped by an inward-migrating planet can eventually be released from resonance 
on polar and even retrograde orbits; this suggests the possibility of non-coplanarity between 
giant planets and much smaller bodies. We think it is likely that the effect seen by the 
authors is actually the M1/M2 limit of the mechanism describe herein (see Section 6). 

Resonant inclination excitation is quite robust against external perturbations. Adding 
the potential of a fairly massive (5 times minimum mass) outer gas disk into the simulation 
does not prevent the onset of the inclination resonance. Neither, in most cases, does the 
addition of a massive inner planetary system. However, it appears that the inclination 
resonances cannot be encountered until the eccentricity of the inner planet gets quite high, 
ei > 0.6. Thus, damping of the eccentricity can prevent excitation of the inclinations. Given 
equal-mass bodies in a 2:1 resonance, we find the threshold relative damping rate for the 
outer body, 62/62, is between 2 and 5 times the relative applied migration rate, d^^^/a2 
(the inner planet is assumed not to interact with the disk). Also, if the mass ratio M-1/M2 
is larger than about 2, the inner body undergoes very little eccentricity excitation and the 
inclination resonances are avoided. 

There is no obvious reason why having two adjacent companions with M1/M2 < 2 should 
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be an uncommon feature of planetary systems. In fact, the majority of discovered multiplanet 
systems, and both of the systems in 2:1 resonance (see below), fulfill this condition at least 
with their projected masses, Msini. Our own Solar System, it should be pointed out, does 
not (Mjupitcr ~ 3Msaturn)- Lcss clcar, though, is the issue of the disk's cfi^cct on planet 
eccentricities. There is not even a consensus on whether the disk-induced eccentricity rate 
of change, Cdisk, is positive or negative. Goldreich & Sari (2003) analytically show that 
the former ought to be the case. However, numerical simulations have tended to support 
the latter; for example. Nelson et al. (2000) perform simulations of gap-opening planets 
with three different hydrodynamic codes and find edisk to be negative in each case. Also, 
in the simulations of Papaloizou, Nelson, & Masset (2001), disk interactions do not cause 
eccentricity growth until the planetary mass exceeds ~ 20 Mj. The problem is that the 
competing effects of eccentricity damping at corotation resonances and eccentricity pumping 
at Lindblad resonances are very nearly of equal strength. Whether one or the other dominates 
depends sensitively on the details of the disk's surface density profile in the vicinity of 
the planet, which is a difficult issue to handle either numerically or analytically. In any 
case, an eccentricity damping timescale of the same order as the orbital decay timescale or 
longer, which we find is required to permit capture into inclination resonances, is entirely 
plausible. A related question is whether the disk might impose some maximum eccentricity 
on embedded planets. As previously mentioned, the standoff distance between a disk and a 
gap-opening planet's orbit is between 3 and 4 Hill radii in numerical simulations (e.g., Bryden 
et al. 2000). Assuming this distance is independent of planetary eccentricity, a Jupiter-mass 
planet will run into the inner disk edge at apocenter when its eccentricity is ~ 0.2 — 0.3, 
which may lead to increased damping. In all of the simulations performed here (except for 
the one shown in Fig. 4, but there the inclination resonance was avoided anyway), the outer 
planet's eccentricity is ~ 0.3 or less at the time the inclination resonance is entered. Thus, 
if such a gap-edge cutoff exists, a pair of resonant planets still has a good chance of entering 
the inclination resonance before suffering significant eccentricity damping. 

We find that once capture into the inclination resonances sets in, the inclinations of the 
planets can rapidly grow to ~ 15° relative to their plane of origin. Since the lines of nodes of 
the two planets' orbits are anti-aligned after entering the inclination resonance, the relative 
inclination is the sum of each orbit's inclination, and grows to ~ 30°. Later, when the system 
undergoes large coupled oscillations in eccentricity and inclination, mutual inclinations can 
get as high as ~ 60° (Fig. 1). Thus, a planetary system which enters inclination resonance 
during its evolution may end up in a highly non-coplanar configuration. This is an important 
possiblity to consider in dynamical investigations of exoplanetary systems; in fact, even the 
inclination relative to our line of sight, and thus the sin i scale factor for the planetary masses, 
may differ significantly among individual members of multi-planet systems. 
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Among the discovered exoplanetary systems, the best candidates for having mutually- 
inchned orbits due to the mechanism described here are those in which the companions are 
(thought to be) in a mean-motion resonance. In the case of the GJ 876 system, Rivera & 
Lissauer (2001) perform dynamical fits allowing for the possiblity of mutual inclinations. 
Most of the systems they find which both provide good fits to the radial velocity data and 
long-term stability have mutual inclinations of 12° or less. They also find one (not necessarily 
the only) stable, locally best-fit (in the sense of a local minimum) region at a very high 
mutual inclination of 77°. Leaving aside this extreme case, since none of the runs performed 
here result in such a high mutual inclination for any significant length of time, one might thus 
conclude that GJ 876 could have undergone resonant inclination excitation, provided some 
mechanism — perhaps eccentricity damping (sec Fig. 6), or direct damping of inclinations 
by the disk, which we have not considered — kept inclinations from getting very high. Also, 
the outer companion's Msini is about 3.4 times that of the inner, so if the outer planet 
is oriented edge-on to our line of sight (sec below), the condition M1/M2 < 2 is fulfilled 
unless the inclination between the two orbits is > 80°. However, the system's relatively low 
eccentricities, < 0.3 for the inner planet, and < 0.1 for the outer (Marcy et al. 2001; Laughlin 
& Chambers 2001), argue against resonant inclination excitation. Unless the eccentricities 
were significantly higher than this in the past, it is unlikely that the inclination resonance 
would ever have been reached. 

The two companions of HD 82943 are also thought to be in a 2:1 resonance, and these 
have inferred eccentricities of 0.41 (outer) and 0.54 (inner) (as announced at an ESO press 
release, April 4, 2001; see http://obswww.unige.ch/ udry/planet/hd82943syst.html). Also, 
the outer companion's M sin i is more than twice that of the inner. This pair of planets may 
therefore provide a more likely setting for resonant inclination excitation. 

The actual mutual inclinations in multiple-planet systems should eventually be revealed 
by astrometric measurements. In fact, the first and thus far only astrometrically-measured 
orbit is that of the outer companion of GJ 876 (Benedict et al. 2002), showing its orbit 
to be oriented nearly edge-on to the line of sight, and thus determining its mass to be 
about 1.9 Mj. The smaller semimajor axis and the (very likely) smaller mass of the inner 
companion, GJ 876c, will make it more difficult to detect astrometrically; if coplanar with 
GJ 876b, it should produce a perturbation in the star's position around 0.2 that produced 
by the outer, or about 0.05 milliarcseconds. Such astrometric precision is attainable in the 
forseeable future; certainly, it lies well within the planned microarcsecond precision of the 
Space Interferometry Mission (SIM). 

A mechanism for tilting planetary orbits relative to the orientation of the original proto- 
planetary disk is particularly interesting in the context of the P Pictoris dust disk (Kalas & 
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Jewitt 1995). The disk, which is oriented nearly edge-on to the hne of sight, displays a warp 
out to a distance of about 50 AU; the disk inside this radius is inclined about 3° relative 
to the outer disk. Mouillet et al. (1997) modeled the warp as arising from the influence of 
a single unseen planet, with a mass between 10~^ and 10^^ that of the central star, with 
a semimajor axis between 1 and 20 AU, and inclined 3 to 5° to the midplane of the outer 
disk. As long as the inclination resonance is actually reached, the simulations we perform 
here indicate that such inclinations could easily arise. Additional structure revealed in re- 
cent observations of P Pic seems even more suggestive: Wahhaj et al. (2003) conclude from 
their observations that the disk in fact contains multiple warps interior to 100 AU, being 
made up of four inchned rings which they denote A, B, C and D, centered at roughly 14, 
28, 52 and 82 AU respectively. They estimate that these these rings are inclined at about 
—32°, +25°, —2° and +2° respectively, relative to the outer disk. They point out that the 
ring radii suggest low-order mean-motion commensurabilities: 3:1 between A and B, and 
2:1 between C and D. There is also the possiblility of a higher-order commcnsurability, 7:3, 
between B and C. Wahhaj et al. (2003) suggest that these rings are associated with planets. 
Their observations appear to fit well with a scenario like the one we explore here, in which 
planets acquire large inclinations relative to the original plane of the protoplanetary disk by 
migrating in resonance with each other. Particularly intriguing are the alternating positive- 
negative inclinations of adjacent rings, mimicing the nodally anti-aligncd configurations that 
the planets take on in our simulations once they acquire significant inclinations. 

Thus far, we have not considered the ultimate fate of our inward-migrating planets. The 
simulations are simply stopped once either planet crosses the inner simulation boundary (0.5 
or 0.1 AU). However, the configuration we are ultimately left with is sensitively dependent 
on how the migration "end game" actually plays out. If migration keeps going down to very 
near the star, one or both planets will be lost. Alternatively, it is possible that the gas disk 
could dissipate sometime after resonant capture but before the planets arc driven into the 
star, though this would require a coincidence in timing. In the absence of other effects that 
outlive the gas, such as interaction with a massive debris disk, interaction with other planets, 
and for bodies with very small periastron distances, tidal interaction with the primary, the 
mean inclinations (and eccentricities) should remain at their values at the time migration 
ceases. Fig. 10 shows an example in which we (abruptly) turn off the induced migration 
of the outer planet at 6 x 10^ years, shortly after the system has entered the inclination 
resonances. 

The longer the migration timescale (i.e., the lower the disk viscosity), the more plausible 
such a coincidence becomes. Another possibility to consider is that migration may under 
some circumstances be outward, for instance if a resonant pair of planets docs not form a 
clean gap and mass flows inward across the gap, as demonstrated by Masset & Snellgrove 
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Fig. 10. — Another version of the basehne run (cf. Fig. 1), this time with migration turned 
off at 5 X 10^ years. 
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(2001). In this case, planets could keep migrating indefinitely without meeting a catastrophic 
end, though they would eventually move beyond the stellocentric radius accessible to radial 
velocity searches. 

Finally, there is the intriguing possibility that inclination excitation could actually serve 
as a "parking mechanism" which puts an end to migration. One should keep in mind that we 
have used a simple, ad-hoc migration prescription in our simulations, nothing more than a 
constant, negative a which is enforced for the outer planet. The interaction of a significantly 
eccentric and inclined giant planet with a disk is in fact highly uncertain. One may speculate 
that a highly- inclined planet might be coupled less strongly to the surrounding disk. At the 
very least, by the time the planet's inclination exceeds that of the disk (~ 4° from midplane 
to scale height at 5 AU, decreasing inward for a flared disk), it may no longer be able to 
maintain a clean gap or edge. Such questions must ultimately be addressed through detailed 
hydrodynamic simulations, which as yet are restricted to the case of planets on circular, 
disk-coplanar orbits. 
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